Each section of this lab has two parts– a model building exercise and a model coding exercise. The material covered here is important and broadly useful – building multi-levels models is a true workhorse for understanding ecological processes because so many problems contain information at nested spatial scales, levels of organization, or categories. It will be worthwhile to dig in deeply to understand it. The big picture is to demonstrate the flexibility that you gain as a modeler by understanding basic principles of Bayesian analysis. To accomplish that, these exercises will reinforce the following:
Ecological data are often collected at multiple scales or levels of organization in nested designs. Group is a catchall term for the upper level in many different types of nested hierarchies. Groups could logically be composed of populations, locations, species, treatments, life stages, and individual studies, or really, any sensible category. We have measurements within groups on individual organisms, plots, species, time periods, and so on. We may also have measurements on the groups themselves, that is, covariates that apply at the upper level of organization or spatial scale or the category that contains the measurements. Multilevel models represent the way that a quantity of interest responds to the combined influence of observations taken at the group level and within the group.
Nitrous oxide N2O, a greenhouse gas roughly 300 times more potent than carbon dioxide in forcing atmospheric warming, is emitted when synthetic nitrogenous fertilizers are added to soils. Qian and colleagues (2010) conducted a Bayesian meta-analysis of N2O emissions (g N \(\cdot\) ha-1 \(\cdot\) d-1) from agricultural soils using data from a study conducted by Carey (2007), who reviewed 164 relevant studies. Studies occurred at different locations, forming a group-level hierarchy (we will use only sites that have both nitrogen and carbon data, which reduces the number of sites to 107 in the analysis here). Soil carbon content (g \(\cdot\) organic C \(\cdot\) g-1 soil dry matter) was measured as a group-level covariate and is assumed to be measured without error. Observations of N2O emission are also assumed to be measured without error and were paired with measurements of fertilizer addition (kg N\(\cdot\) ha-1 \(\cdot\) year-1). The effect of different types of fertilizer was also studied.
You are going to use these data to build increasingly complex models of N2O emission. The initial models will ignore some important covariates as well as how the data are structured hierarchically into sites. This is ok! When writing for a multi-level model like this one, do it incrementally, starting with a separate model for each site (the no-pool model) or a model that ignores sites entirely (the pooled model). After getting these models to work you can add complexity by drawing the intercept for each model from a distribution, before pursuing further refinements. We strongly sugggest this approach because it is always best to do the simple thing first: there is less to go wrong. Also, when things do go wrong it will be clearer as to what is causing the problem.
You need to load the following libraries. Set the seed to 10 to compare your answers to ours.
library(actuar)
library(rjags)
library(ggplot2)
library(ggthemes)
library(gridExtra)
library(MCMCvis)
library(HDInterval)
library(SESYNCBayes)
library(plyr)
library(tidyverse)
set.seed(10)
Let’s begin by ignoring the data on soil carbon, site, and fertilizer type so that all observations are drawn from a single pool. This is what’s known as complete pooling (see Gelman and Hill, 2007), or just a pooled model. You will use a linearized power function for your deterministic model of emissions as a function of nitrogen input:
\[ \begin{aligned} \mu_{i} & = \gamma x_{i}^{\beta}\\ \alpha & = \log\big(\gamma\big)\\ \log\big(\mu_{i}\big) & = \alpha+\beta\big(\log(x_i)\big)\\ g\big(\alpha,\beta,\log(x_i)\big) & = \alpha+\beta\big(\log(x_i)\big) \\ \end{aligned} \]
It is always a good idea to look at the data. Examine the head of the data frame for emissions. Note that the columns group.index and fert.index contain indices for sites and fertilizer types. We are going to ignore these for now since the pooled model does not take these into account. Use the code below to plot N2O emissions as a function of fertilizer input for both the logged and unlogged data.
head(N2OEmission)
## fertilizer group carbon n.input emission reps group.index fert.index
## 1 A 14 2.7 180 0.620 13 10 2
## 2 A 14 4.6 180 0.450 13 10 2
## 3 A 11 0.9 112 0.230 12 7 2
## 4 A 38 0.5 100 0.153 14 29 2
## 5 A 1 4.0 250 1.000 6 1 2
## 6 A 38 0.5 100 0.216 14 29 2
We are going to use ggplot to visualize the data in this lab. If you are unfamiliar with this package, don’t worry. We will provide you will all the codes you need and help your get oriented. We think you will find the plotting functions in ggplot very powerful and intuitive. We start by using ggplot to load the data frame we will plot data from. Then we add geom_point and use the aes argument (the aesthetic mappings) to define the x and y values for the points. All ggplot functions require you to define the aesthetic mappings as needed. Here, they are the same as setting x and y in the normal plot functions. The other big difference is that ggplot allows you to add successive layers to the plot using the + operator. You will see later on that this offers a lot of flexibility. We add the geom_line feature and then set the theme to minimal. Lastly, we use the grid.arrange function to position multiple plots at once. This is similar to using mfrow with par.
g1 <- ggplot(data = N2OEmission) +
geom_point(aes(y = emission, x = n.input), alpha = 3/10, shape = 21, colour = "black",
fill = "brown", size = 3) +
theme_minimal()
g2 <- ggplot(data = N2OEmission) +
geom_point(aes(y = log(emission), x = log(n.input)), alpha = 3/10, shape = 21, colour = "black",
fill = "brown", size = 3) +
theme_minimal()
gridExtra::grid.arrange(g1, g2, nrow = 1)
You will now write a simple, pooled model where you gloss over differences in sites and fertilizer types and lump everything into a set of \(x\) and \(y\) pairs using the R template provided below. It is imperative that you study the data statement and match the variable names in your JAGS code to the left hand side of the = in the data list. Call the intercept alpha, the slope beta and use sigma to name the standard deviation in the likelihood. Also notice, that we center the nitrogen input covariate to speed convergence. You could also standardize this as well.
In addition to fitting this model, we would like you to have JAGS predict the mean logged N2O emissions and the median unlogged N2O emissions as a function of soil fertilizer input. (Why median? Hint: think back to the distribution of the untransformed data above in question 3 above). To help you out we have provided the range of N2O values to predict over as the third element in the data list. Make sure you understand how we chose these values.
Note that in this problem and the ones that follow we have set up the data and the initial conditions for you. This will save time and frustration, allowing you to concentrate on writing code for the model but you must pay attention to the names we give in the data and inits lists. These must agree with the variable names in your model. Please see any of the course instructors if there is anything that you don’t understand about these lists.
n.input.pred <- seq(min(N2OEmission$n.input), max(N2OEmission$n.input), 10)
data = list(
log.emission = log(N2OEmission$emission),
log.n.input.centered = log(N2OEmission$n.input) - mean(log(N2OEmission$n.input)),
log.n.input.centered.pred = log(n.input.pred) - mean(log(N2OEmission$n.input)))
inits = list(
list(alpha = 0, beta = .5, sigma = 50),
list(alpha = 1, beta = 1.5, sigma = 10),
list(alpha = 2, beta = .75, sigma = 20))
Let’s overlay the predicted mean logged N2O emissions and median unlogged N2O emissions as a function of soil fertilizer input from the pooled model on top of the raw data. We summarize the predictions using MCMCpstr() twice - once to get the 95% HDPI intervals and a second time to get the posterior median for each fertilizer input value. We combine these predictions into two data frames, one for the logged N2O emissions and one for untransformed N2O emissions. We append our new graphical elements onto our old plots with the + operator. We plot the median of the posterior distribution as a black line with geom_line() and the 95% credible intervals as a yellow shaded region using the geom_ribbon() function. These data come from a different data frame than the one we used to plot the raw data so we need to add the data argument in the new geom_line and geom_ribbon. Again, we provide you with the code to do this to save time. You will need to modify this code to make similar plots for models you fit in later exercises.
pred_po1 <- MCMCpstr(zc.pooled, params = c("mu_pred", "log_mu_pred"), func = function(x) hdi(x, .95))
pred_po2 <- MCMCpstr(zc.pooled, params = c("mu_pred", "log_mu_pred"), func = median)
pred.po.df <- cbind(n.input.pred, data.frame(pred_po1$mu_pred), mid = pred_po2$mu_pred)
lpred.po.df <- cbind(log.n.input.pred = log(n.input.pred), data.frame(pred_po1$log_mu_pred), mid = pred_po2$log_mu_pred)
g3 <- g1 +
geom_line(data = pred.po.df, aes(x = n.input.pred, y = mid)) +
geom_ribbon(data = pred.po.df, aes(x = n.input.pred, ymin = lower, ymax = upper), alpha = 0.2, fill = "yellow")
g4 <- g2 +
geom_line(data = lpred.po.df, aes(x = log.n.input.pred, y = mid)) +
geom_ribbon(data = lpred.po.df, aes(x = log.n.input.pred, ymin = lower, ymax = upper), alpha = 0.2, fill = "yellow")
gridExtra::grid.arrange(g3, g4, nrow = 1)
Great! - you’ve got the pooled model fitted and made some predictions from it. However, perhaps the idea of ignoring the site effects is not sitting so well with you. Let’s take this a step further by modeling the relationship between N2O emission and fertilizer input such that the intercept \(\alpha\) varies by site (we will again ignore the data on soil carbon and fertilizer type). This is the opposite of the pooled model where we completely ignored the effect of site. Here we treat the intercept for each site as independent. This is commonly called a no-pool model.
\[ \begin{aligned} \mu_{ij} & = \gamma x_{ij}^{\beta}\\ \alpha_{j} & = \log\big(\gamma_{j}\big)\\ \log\big(\mu_{ij}\big) & = \alpha_{i}+\beta\big(\log(x_{ij})\big)\\ g\big(\alpha_{j},\beta,\log(x_{ij})\big) & = \alpha_{j}+\beta\big(\log(x_{ij})\big) \\ \end{aligned} \]
Let’s visualize the data again, but this time highlighting the role site plays in determining the relationship between N2O emission and fertilizer input. First, head() the data to see how groups are organized. You will use group.index to group the observations by site.
head(N2OEmission)
## fertilizer group carbon n.input emission reps group.index fert.index
## 1 A 14 2.7 180 0.620 13 10 2
## 2 A 14 4.6 180 0.450 13 10 2
## 3 A 11 0.9 112 0.230 12 7 2
## 4 A 38 0.5 100 0.153 14 29 2
## 5 A 1 4.0 250 1.000 6 1 2
## 6 A 38 0.5 100 0.216 14 29 2
Use the code below to plot logged N2O emissions against logged fertilizer input. This is the same ggplot code as before except now we amend it to make plots for individual sites simply by adding the facet_wrap function and specifying the grouping variable(here it is group.index) as an argument.
g2 + facet_wrap(~group.index)
You will now write a simple, no-pool model using the R template provided below. Again notice, that we center the nitrogen input covariate to speed convergence. You could also standardize this as well.
In addition to fitting this model, we would like you to have JAGS predict the mean logged N2O emissions for each site as a function of soil fertilizer input. To help you out we have provided the range of N2O values to predict over as the third element in the data list.
n.sites <- length(unique(N2OEmission$group.index))
n.input.pred <- seq(min(N2OEmission$n.input), max(N2OEmission$n.input), 10)
data = list(
log.emission = log(N2OEmission$emission),
log.n.input.centered = log(N2OEmission$n.input) - mean(log(N2OEmission$n.input)),
log.n.input.centered.pred = log(n.input.pred) - mean(log(N2OEmission$n.input)),
group = N2OEmission$group.index,
n.sites = n.sites)
inits = list(
list(alpha = rep(0, n.sites), beta = .5, sigma = 50),
list(alpha = rep(1, n.sites), beta = 1.5, sigma = 10),
list(alpha = rep(-1, n.sites), beta = .75, sigma = 20))
We modify the MCMCpstr code from the previous model to produce a data frame of the median and 95% HDPI credible intervals of N2O emission predictions for each site. MCMCpstr preserves the shape of the parameter from your JAGS model, which can be very handy in certain situations. Here, pred_np1 is a 3D-array, whose rows are fertilizer inputs, columns are sites, and whose z-values are the quantities produced by the hdi function, which in this case is the lower and upper credible interval. You can str the pred_np1 object to see this for yourself. For plotting purposes though, we would like a data frame with columns for site, fertilizer input, the posterior’s median emission, and the posterior’s lower and upper HDPI credible intervals. This can be made easily using the adply function. We use cbind to make the data frame we seek, taking advantage of the fact that n.input.pred will repeat each site, which is exactly what we want it to do. Lastly, we rename X2 (created during the adply step) to group.index for clarity.
pred_np1 <- MCMCpstr(zc.nopooled, params = "log_mu_site_pred", func = function(x) hdi(x, .95))
pred_np2 <- MCMCpstr(zc.nopooled, params = "log_mu_site_pred", func = median)
pred_np1_df <- plyr::adply(pred_np1$log_mu_site_pred, c(1, 2))[, 2:4]
pred_np2_df <- plyr::adply(pred_np2$log_mu_site_pred, c(1, 2))[, 3]
lpred.np.df <- cbind(log.n.input.pred = log(n.input.pred), pred_np1_df, median = pred_np2_df) %>%
rename(group.index = X2)
To add the predictions to the plots for each site we use geom_line and geom_ribbon again, in combination with facet_wrap.
g2 +
geom_line(data = lpred.np.df, aes(x = log.n.input.pred, y = median)) +
geom_ribbon(data = lpred.np.df, aes(x = log.n.input.pred, ymin = lower, ymax = upper), alpha = 0.2, fill = "yellow") +
facet_wrap(~group.index)
So far you have either ignored sites completely (the pooled model) or treated all the sites intercepts as independent from one another (the no-pool model). Now you are going to treat the sites as partially pooled, meaning you will model the site intercepts as coming from a common distribution. In other words, you will treat the intercept in your model as a group level effect (aka, random effect). The deterministic model of N2O emissions remains a linearized power function, but two subscripts are required: \(i\) indexes the measurement within sites and \(j\) indexes site. Assume that the intercepts are drawn from a distribution with mean \(\mu_{\alpha}\) and variance \(\varsigma_{\alpha}^2\).
Now you will implement the model that allows intercept to vary by group, where each intercept is drawn from a common distribution. Again, use the template provided below to allow you to concentrate on writing JAGS code for the model. Note that you must use the index trick covered in lecture to align the different groups with different intercepts. Here are the preliminaries to set up the model:
n.input.pred <- seq(min(N2OEmission$n.input), max(N2OEmission$n.input), 10)
n.sites <- length(unique(N2OEmission$group.index))
data = list(
log.emission = log(N2OEmission$emission),
log.n.input.centered = log(N2OEmission$n.input) - mean(log(N2OEmission$n.input)),
log.n.input.centered.pred = log(n.input.pred) - mean(log(N2OEmission$n.input)),
group = N2OEmission$group.index,
n.sites = n.sites)
inits = list(
list(alpha = rep(0, n.sites), beta = .5, sigma = 50, mu.alpha= 0, sigma.alpha = 10),
list(alpha = rep(1, n.sites), beta = 1.5, sigma = 10, mu.alpha= 2, sigma.alpha = 20),
list(alpha = rep(-1, n.sites), beta = .75, sigma = 20, mu.alpha= -1, sigma.alpha = 12))
Use the code from the pooled model to visualize the model predictions again.
In the previous example, we assumed that the variation in the intercept was attributable to spatial variation among sites. We did not try to explain that variation, we simply acknowledged that it exists. Now we are going to “model a parameter” using soil carbon content data at the site-level to explain variation in the intercepts among sites. Modify the previous model to represent the effect of soil carbon on the intercept using the deterministic model. Here, we logit transform the carbon data to “spread them out” mapping 0-1 to all real numbers.
\[g_2\big(\kappa,\eta,\text{logit}(w_j)\big) =\kappa + \eta \text{logit}w_j\] to predict \(\alpha_j\).
Modify your random intercepts model to implement the model that include soil carbon content as covariate at the site level. Make predictions for how mean logged N2O emission and median N2O emission varies with respect to soil fertilizer input across all possible sites of average soil carbon content. We again provide you with the data and initial values.
n.input.pred <- seq(min(N2OEmission$n.input), max(N2OEmission$n.input), 10)
n.sites <- length(unique(N2OEmission$group.index))
data = list(
log.emission = log(N2OEmission$emission),
log.n.input.centered = log(N2OEmission$n.input) - mean(log(N2OEmission$n.input)),
log.n.input.centered.pred = log(n.input.pred) - mean(log(N2OEmission$n.input)),
w = log(SiteCarbon$mean / (100 - SiteCarbon$mean)),
group = N2OEmission$group.index,
n.sites = n.sites)
inits = list(
list(alpha = rep(0, n.sites), beta = .5, sigma = 50, sigma.alpha = 10, eta = .2, kappa = .5),
list(alpha = rep(1, n.sites), beta = 1.5, sigma = 10, sigma.alpha = 20, eta = 3, kappa = .7),
list(alpha = rep(-1, n.sites), beta = .75, sigma = 20, sigma.alpha = 12, eta = .1, kappa = .3))
Use the code from the pooled and no-pooled models to visualize the model predictions again. How does modeling site soil carbon affect the uncertainty in predicting N2O emissions for new sites?